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ABSTRACT 

We analyze the mass accretion histories (MAHs) and density profiles of cluster-size halos with virial masses 
of 0.6-2.5 X 10'"^/!"' Mq in a flat ACDM cosmology. We find that most MAHs have a similar shape: an early, 
merger-dominated mass increase followed by a more gradual, accretion-dominated growth. For some clusters the 
intense merger activity and rapid mass growth continue until the present-day epoch. In agreement with previous 
studies, we find that the concentration of the density distribution is tightly correlated with the halo's MAH and with 
its formation redshift. During the period of fast mass growth the concentration remains approximately constant 
and low Cy « 3 - 4, while during the slow accretion stages the concentration increases with decreasing redshift 
as Cv oc (1 +z)~'. We consider fits of three widely discussed analytic density profiles to the simulated clusters 
focusing on the most relaxed inner regions. We find that there is no unique best fit analytic profile for all the 
systems. At the same time, if a cluster is best fit by a particular analytic profile at z = 0, the same is usually true at 
earlier epochs out to z ^ 1 -2. The local logarithmic slope of the density profiles at 3% of the virial radius ranges 
from -1 .2 to -2.0, a remarkable diversity for the relatively narrow mass range of our cluster sample. Interestingly, 
for all the studied clusters the logarithmic slope becomes shallower with decreasing radius without reaching an 
asymptotic value down to the smallest resolved scale (< 1 % of the virial radius). We do not find a clear correlation 
of the inner slope with the formation redshift or the shape of the halo's MAH. We do find, however, that during 
the period of rapid mass growth the density profiles can be well described by a single power law p(r) oc with 
7 ~ 1.5-2. The relatively shallow power law slopes result in low concentrations at these stages of evolution, 
as the scale radius where the density profiles reaches the slope of -2 is at large radii. This indicates that the 
inner power law like density distribution of halos is built up during the periods of rapid mass accretion and active 
merging, while outer steeper profile is formed when the mass accretion slows down. To check the convergence 
and robustness of our conclusions, we resimulate one of our clusters using eight times more particles and twice 
better force resolution. We find good agreement between the two simulations in all of the results discussed in our 
study. 

Subject headings: cosmology: theory - dark matter - clusters: formation - clusters - structure methods: 
numerical 



1. INTRODUCTION 

During the last decade, there has been an increasingly grow- 
ing interest in testing the predictions of variants of the cold dark 
matter models (CDM) on small scales. The interest was spurred 
by indications that the density distribution in the inner regions 
of dark matter halos predicted by CDM is at odds with the ob- 
served galactic rotation curves (Flores & Primack 1994; Moore 
1994). This discrepancy is yet to be convincingly resolved and 
is still a subject of active debate (e.g.. Cote et al. 2000; van den 
Bosch & S waters 2001; Blais-Ouellette et al. 2001; de Blok 
et al. 2001, 2003; Swaters et al. 2003). In addition, the CDM 
models face other apparent discrepancies with observations on 
galactic scales such as the amount of substructure in galactic 
halos (Klypin et al. 1999; Moore et al. 1999a), the incorrect 
normalization of the TuUy-Fisher relation, the angular momen- 
tum of disk galaxies (Navarro & Steinmetz 1997, 2000), the 
ellipticity of dark matter halos (Ibata et al. 2001), and others. 

In the past several years, the density distribution in the cores 
of galaxy clusters has also become a subject of a related de- 
bate. CDM models predict cuspy density profiles without flat 
cores (Frenk et al. 1985; Quinn et al. 1986; Dubinski & Carl- 
berg 1991). NavaiTO et al. (1996, 1997, NEW) argued that the 
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CDM halo profiles can be described by the following simple 
formula in all cosmologies and at all epochs, 

Po 



pir): 



(1) 



(r/r,)(l + r/r,)2' 
This analytic formula describes the density profile of a halo us- 
ing two parameters: a characteristic density, po, and a scale 
radius, r,. These parameters are determined by the halo virial 
mass. My, and concentration index, c = r^/rs, where r,, is the 
virial radius of the halo. In addition, NEW argued that there is 
a tight correlation between c and My, which means that the halo 
profiles of different mass objects form a one parameter family. 

Moore et al. (1998) (see also Ghigna et al. 2000) carried out 
a convergence study of the dark matter profiles and concluded 
that high mass resolution is required to resolve the inner den- 
sity distribution robustly. They advocated the analytic density 
profile of the form p{r) oc (r/rj)"''^[l + (r/r,)''^]"^ as a bet- 
ter description of the density distribution of their simulated ha- 
los. This profile behaves similarly to the NEW profile at large 
radii (oc r~^), but is steeper at small radii (oc r"'-^). Fukushige 
& Makino (1997, 2001, 2003) reached similar conclusions us- 
ing a set of independent simulations. Jing & Suto (2000) pre- 
sented a systematic study of the density profiles of halos with 
masses in the range 2 x 10'^ -5 x lO'"^ Mq. They found 
that the inner slope at a radius of 1% of the virial radius is shal- 
lower (« -1.1) for cluster mass halos than for galactic halos 
(« -1 .5). Recently, Hayashi et al. (2003); Navarro et al. (2003) 
found that often the logarithmic slope of the density distribu- 
tion at the convergence radius is steeper than -1 as expected 
from the NEW profile, but significantly shallower than the -1.5 
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inner slope found by Moore et al. (1998). Several other studies 
(Kravtsov et al. 1997, 1998; Avila-Reese et al. 1999; Jing 2000; 
Bullock et al. 2001; Klypin et al. 2001; Fukushige et al. 2003) 
found a significant scatter in both the shape of the density pro- 
files and halo concentrations, likely related to the details of the 
mass accretion histories of individual objects (Wechsler et al. 
2002; Zhao et al. 2003b). Mucket & Hoeft (2003) and Hoeft 
et al. (2004) study the radial dependence of the gravitational 
potential and the velocity dispersion and come to the conclu- 
sion that there does not exist a slope asymptote of the density 
profile over a wide range but the slope increases with decreas- 
ing radius and reaches the value -0.58 as r ^ 0. 

Observational constraints on the inner slope of the dark mat- 
ter density distribution in galactic halos are difficult because the 
distribution is affected by the cooling and dynamics of baryons. 
The dark matter profiles in clusters, on the other hand, should 
be less affected by cooling as a much smaller fraction of cluster 
baryons is observed to be in the cold condensed phase. Ob- 
servational studies of the mass distribution in clusters using 
weak lensing and hydrostatic analysis of the X-ray emitting 
gas show that the overall mass distribution is in general agree- 
ment with CDM predictions (Allen 1998; Clowe et al. 2000; 
Willick & Padmanabhan 2000; Clowe & Schneider 2001 ; Shel- 
don et al. 200 1 ; Arabadjis et al. 2002; Athreya et al. 2003 ; Bautz 
& Arabadjis 2003). 

Strong lensing studies can probe the mass distribution in the 
inner region of clusters and thus test the "cuspiness" of clus- 
ter halos. However, the results of strong lensing analyses have 
so far been contradictory, even in the case when the same sys- 
tem was studied. Tyson et al. (1998), for example, argue that 
the density profile of cluster CL0024+1654 has a constant den- 
sity core, while Broadhurst et al. (2000) find that the mass dis- 
tribution in this cluster is cuspy. Czoske et al. (2002) argue 
that CL0024+1654 is undergoing a major merger and its den- 
sity profile may not be representative. Sand et al. (2002) find 
that the inner slope of the density profile in cluster MS2137- 
23 is flatter than expected in CDM models, a conclusion they 
recently confirmed for six more clusters (Sand et al. 2003). 
Gavazzi et al. (2002), reanalyzing the same observations, ar- 
gue that if the fifth demagnified image near the center of the 
lensing potential is not taken into account then the inner slope 
may be consistent with CDM predictions. 

Given the disagreement among the different analytical fits 
proposed for the density profiles of dark matter halos found in 
simulations and a possible discrepancy with strong lensing ob- 
servations, it is interesting to conduct a systematic study of the 
density profiles of clusters in the concordance ACDM model. 
The study of cluster mass halos is also interesting because the 
typical concentrations of their matter distribution are lower than 
those of galactic halos. Thus, if an asymptotic inner slope, sug- 
gested by the analytic profiles, does exist it should be reached at 
a larger fraction of the virial radius in cluster halos and should 
be easier to detect. 

The paper is organized as follows. In the following two sec- 
tions, we describe the numerical simulations and halo finding 
algorithm used in our analysis. In § 4 we discuss the mass ac- 
cretion histories of the analyzed clusters. In § 5 we present the 
convergence test, discuss the fitting procedure, and our results 
on the shapes and inner slopes of the density profiles. We sum- 
marize our results and conclusions in § 6. 

2. NUMERICAL SIMULATIONS 



We use the Adaptive Refinement Tree code (Kravtsov et al. 
1997) to follow the evolution of cluster-size halos in the flat 
ACDM cosmology: (Om,^^A,/j,cr8)= (0.3,0.7,0.7,0.9). We use 
the initial spectrum in the Holtzman approximation with flb = 
0.03 (see Klypin & Holtzman 1997). The code starts with a 
uniform 256^ grid covering the entire computational box. This 
grid defines the lowest (zeroth) level of resolution. Higher force 
resolution is achieved in the regions corresponding to collaps- 
ing structures by recursive adaptive refinement of all such re- 
gions. Each cell can be refined or de-refined individually. The 
cells are refined if the particle mass contained within them ex- 
ceeds a certain specified threshold value. The code thus refines 
to follow the collapsing objects in a quasi-lagrangian fashion. 

The cluster halos were simulated in a box of 80/z~' Mpc. A 
low resolution simulation was run first. A dozen cluster ha- 
los were identified and multiple mass resolution technique was 
used to set up initial conditions (Klypin et al. 2001). Namely, a 
lagrangian region corresponding to a sphere of radius equal to 
two virial radii around each halo was re-sampled with the high- 
est resolution particles of mass Wp = 3.16 x lO^^^'M©, corre- 
sponding to an effective number of 5 12-' particles in the box, at 
the initial redshift of the simulation (zi = 50). The high mass res- 
olution region was surrounded by layers of particles of increas- 
ing mass with a total of three particle species. Only regions 
containing highest resolution particles were adaptively refined 
and the threshold for refinement was set to correspond to a mass 
of 4 highest resolution particles per cell. Each cluster halo is re- 
solved with ^ lO*" particles within its virial radius at z = 0. The 
size of the highest refinement level cell was 1 .2h~^ kpc. In addi- 
tion, one of the clusters was re-simulated with eight times more 
particles (wp = 3.95 x 1O^/!~'M0) to study the convergence of 
the density profiles. In this simulation, the smallest cell size 
reached was 0.6h~^ comoving kpc. 

The time steps were chosen so that no particle moves by more 
than a fraction of the parent cell size in a single step. This cri- 
terion was motivated by the convergence studies presented by 
Klypin et al. (2001). For the analyzed simulations, the num- 
ber of steps at the highest refinement level was w 250000 or 
Af « 2-3 X 10"* yrs and a factor of 2 larger for each lower 
refinement level. For the high-resolution resimulation of one 
of the cluster used for convergence check, the number of steps 
at the highest refinement level was « 500000. We analyze the 
cluster profiles and their mass accretion histories using 19 out- 
puts from z= 10 to z = 0, with a typical time interval between 
outputs of ~ 0.7 Gyr. 

3. HALO IDENTIFICATION 

To identify cluster halos we use a variant of the Bound Den- 
sity Maxima (BDM) halo finding algorithm. The main idea of 
the BDM algorithm is to find positions of local maxima in the 
density field smoothed at a certain scale and to apply physically 
motivated criteria to test whether the identified site corresponds 
to a gravitationally bound halo. The detailed description of the 
algorithm is given in Klypin & Holtzman (1997) and Klypin 
et al. (1999). 

We start by calculating the local overdensity at each parti- 
cle position using the SPH smoothing kernel"* of 24 particles. 
We then sort particles according to their overdensity and use all 
particles with 6 > Smin = 5000 as potential halo centers. Starting 
with the highest overdensity particle, we surround each poten- 

To calculate the density we use the pubUcly available code smooth: 
http : //www-hpcc .astro . Washington . edu/tools/tools .html 
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Table 1 
Simulated Halo parameters 



Halo 


Mm 


nso 


Vmax 




Qf' M@) 


(r'Mpc) 


(kms"') 


CLl 


2.5 X lO''* 


1.58 


973 


CL2 


2.4 X lO"* 


1.56 


1011 


CL3 


2.3 X lO'" 


1.55 


904 


CL4 


1.3 X 10"* 


1.29 


826 


CL5 


1.3 X 10" 


1.27 


798 


CL6 


1.2 X 10" 


1.25 


785 


CL7 


1.2 X 10" 


1.23 


587 


CL8 


1.2 X 10" 


1.23 


695 


CL9 


9.7 X 10" 


1.16 


630 


CLIO 


8.6 X 10" 


1.11 


597 


CLU 


8.1 X 10" 


1.09 


670 


CL12 


8.1 X 10" 


1.09 


758 


CL13 


7.3 X 10" 


1.05 


607 


CL14 


5.8 X 10" 


0.98 


603 



tial center by a sphere of radius rf^d = 50h kpc and exclude 
all particles within this sphere from further center search. Af- 
ter all the potential centers are identified, we analyze the den- 
sity distribution and velocities of the surrounding particles to 
test whether the center corresponds to a gravitationally bound 
clump (Klypin et al. 1999). We then construct profiles using 
only bound particles and use them to calculate the properties 
of halos such as the maximum circular velocity V^ax, the mass 
M, etc. In this study, we consider only isolated cluster-size ha- 
los. We should note that for isolated halos the BDM algorithm 
works very similarly to the commonly used spherical overden- 
sity (SO) algorithm. 

The virial radius is a convenient measure of the halo size. We 
define the virial radius as the radius within which the density is 
equal to 180 times the average density of the universe at a given 
epoch. The separation between halos is sometimes smaller than 
the sum of their virial radii. In such cases, the definition of the 
outer boundary of a halo and its mass are somewhat ambigu- 
ous. To this end, in addition to the virial radius, we estimate 
the truncation radius, r^, at which the logarithmic slope of the 
density profile constructed from the bound particles becomes 
larger than -0.5 as we do not expect the density profile of the 
CDM halos to be flatter than this slope. In general we consider 
the halo radius to be = min(ri8o, ''t)- 

In our analysis we use only clusters with masses 
> 5 X lO'-'/z"' Mq. In most cases two or more clusters were 
identified with this mass threshold in each run. To distinguish 
between the isolated cluster halos and massive subhalos we use 
additional information, such as the virial-to-tidal radius ratio, 
the maximum circular velocity, and the number of gravitation- 
ally bound particles within rt,. We consider halos to be isolated 
if their separation is larger than one third of the sum of their 
virial radii. We list the present-day properties of the cluster 
halos included in our sample in Table 1 . The masses and radii 
correspond to the cumulative overdensity of 180 times the mean 
density of the Universe. In this list, clusters 4/12 and 10/13 are 



close pairs, clusters 6/11/14 and 7/8/9 are triplets, while clus- 
ters 1, 2, 3, and 5 are weU-isolated systems. The clusters thus 
sample a variety of environments. 

We stress that the clusters in the analyzed sample were se- 
lected randomly - no specific criterion of relaxation or sub- 
structure was used. As the Figure 1 below shows clusters in 
our sample span a wide range of mass accretion histories and 
formation redshifts. 

4. MASS ACCRETION HISTORIES 

In the hierarchical structure formation scenario halos are as- 
sembled via a continuous process of merging and accretion. 
Details of the mass accretion history (MAH) may affect the 
shape of the halo and its density distribution (Navarro et al. 
1997; Bullock et al. 2001; Wechsler et al. 2002; Zhao et al. 
2003b). It is therefore interesting to study the accretion history 
of the halos in conjunction with the study of their density pro- 
files. In this section we study the details of the assembly of the 
simulated cluster halos by following the most massive progen- 
itor from z = 10 up to the present. We will discuss connections 
between the halo density profile and its MAH in § 5.3 and 6. 

4.1. Constructing MAHs 

For each z = cluster halo we identify the most massive pro- 
genitor using the halo catalogs of the previous time output. In 
what follows, if halo 1 is the most massive progenitor of halo 
2, then halo 2 will be referred to as the offspring of halo 1. To 
identify the most massive progenitor of a halo we first iden- 
tify all of its progenitors in the halo catalog. We then eUminate 
from the set of potential progenitors, objects that are signifi- 
cantly tidally stripped (i.e., their virial-to-tidal radius ratio is 
greater than 3.5). We use the following criteria to identify the 
most massive progenitor among the remaining candidates. 1. 
We eliminate candidate progenitors with masses less than 20% 
of the offspring mass. 2. Using the peculiar velocity of the 
offspring, we find the approximate location of the progenitor 
in the previous output and eliminate the candidates outside the 
sphere with radius r = 10 x VpAt, where Vp is the offspring pe- 
culiar velocity and At the time elapsed between two successive 
outputs. 3. We require that candidate progenitor and offspring 
halos have a certain fraction of common particles. In what fol- 
lows, /i (f2) denotes the ratio of the number of particles that 
offspring and progenitor have in common to the number of par- 
ticles in the offspring (progenitor). 

The candidate with the largest number of common particles 
with the offspring and with /i > 0.5 is then chosen to be the 
most massive progenitor. At the same time the condition /2 > 
0.5 is also checked and found to be satisfied. If all the progeni- 
tors have /i < 0.5, as is often the case during major mergers, for 
the progenitor with the largest /i we also require that more than 
95% of the particles within a comoving radius of 10 kpc 
from the most bound particle of the progenitor are also found 
in the offspring. Starting from z = 0, we repeat the identification 
of the most massive progenitor for all the 19 simulation output 
epochs to z — 10, or until the progenitor can no longer be iden- 
tified. The mass accretion histories constructed in this way for 
each of the analyzed clusters are shown in Figure 1 (solid lines). 

Most of the MAHs have a qualitatively similar shape: a rapid 
increase in mass during the early epochs and a relatively slow 
increase at the later stages of evolution. Despite the similar- 
ities, the details of MAHs differ significantly from object to 
object. CL9, CLIO, and CL13 have not yet reached the second. 
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Fig. 1. — Mass accretion histories of the cluster halos {solid lines). Also shown are the analytic fits of Eq. (4), M(a)/Mo = a''exp[-o:(a- 1)] (where a = a/ao, and 
a = (1 +z)"'). The dashed lines show fits with both a and p varied, while the dotted lines show the fits with the parameter p fixed to zero. The formation redshift 
Zf, given by Eq. (3) in terms of a, is shown in the legend of each panel. The best fit a obtained from the two fits is nearly identical in all cases, except for CL9 and 
CL13. For CL9 and 13 we also plot the fit obtained using Eq. (5) (crosses). Finally, to show the effect of a recent major merger, for CL7 we plot the fit assuming an 
epoch of observation ao = 0.95 rather than oq = 1. used for aU the other fits (dash-dotted line). The Zf value obtained in this case is given in the parentheses. 
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accretion dominated stage of their evolution. Their mass is ac- 
cumulated via intense merger activity up to the present epoch. 
CL7 and CL8 appear to have reached the slow accretion phase, 
but experienced a late major merger. The masses of CLll and 
14 at early epochs increase almost linearly with the expansion 
factor (log (M/Mq) oc a), and reach a M w const plateau at later 
epochs. 

4.2. Major mergers 

In addition to the overall shape of the halo MAH, it is useful 
to have some more specific information on the major mergers 
experienced by the clusters. We will use the term major merger 
to describe all events that result in a more than 30% increase in 
the mass of the main progenitor between the two output epochs. 
As we mentioned above, the average time elapsed between two 
successive outputs of the simulation is w 0.7 Gyr, which is close 
to the crossing time of w 1 Gyr for a wide range of halo masses. 
The crossing time is a lower limit for the merger time scale, 
which means that the spacing of our outputs is appropriate for 
merger identification (see also tests in Gottlober et al. 2001). 
We tabulate the redshift of the last major merger, Zlmm, as well 
as the corresponding fractional mass change, AM/M, for all 
the clusters in columns 2 and 3 of Table 2, respectively. One 
should keep in mind that these numbers are only indicative, 
since defining a major merger, e.g., as a 20% mass increase, 
would render zlmm for CL4 equal to ~ 0.15. 

4.3. Formation redshift and MAH shape 

To characterize evolution of the halos, one can introduce the 
halo formation epoch (or redshift). Usually, the formation epoch 
is defined as the time when the mass in the most massive pro- 
genitor(s) is equal to some fraction of the halo's final mass. Mo 
(e.g., see Lacey & Cole 1993; Navarro et al. 1997). Taking this 
fraction to be equal to 1 /2 we calculate the formation redshift, 
zi/2, which we tabulate in column 4 of Table 2. In order to find 
z 1/2 we use linear interpolation between the successive outputs 
that bracket M/Mo = 1/2. It is interesting to note that Z\i2 is 
typically smaller than zlmm- 

As pointed out by Wechsler et al. (2002, hereafter W02), 
defining the formation redshift as the redshift where the ratio 
M/Mq takes a specific value gives a formation redshift that de- 
pends on the time of observation of the halo. In addition, the 
definition uses the MAH of the halo at two epochs only (the 
formation and present epochs) and therefore makes it sensitive 
to the local jumps in the MAH and less sensitive to the overall 
MAH shape. The case of CL7 may serve as an illustration. This 
object had entered its quiescent stage of evolution relatively 
early. Nevertheless, the low value of its formation redshift, zi/2, 
is determined largely by the single late major merger. In view of 
these considerations, W02 argued that the mass accretion his- 
tories can be better characterized by a formation redshift that is 
derived from a functional fit to the entire MAH. Namely, they 
propose to fit the MAHs of halos by a simple exponential: 

M(fl) = exp [a (l - 1/fl)] ; d = a/ao; a = (l+z)~', (2) 

where M = M/Mq, and Mq and ao are the virial mass of the halo 
and expansion factor at the epoch of observation, respectively. 
Using the fit, one can define the formation epoch independent of 
the epoch of observation as the redshift corresponding to a fixed 
value of d logM/ d loga = S. The value of S is arbitrary, and we 
follow W02 and choose S = 2 since this is the value required to 
match the concentration index-collapse redshift relation found 



by Bullock et al. (2001). The formation redshift, Zf, can then 
be defined by the relation 



z/ = -(1+zo)- 

a 



1. 



(3) 



The formation redshift in this definition is independent of the 

epoch of observation; the factor (1 +zo) appears only in order 
to comply with the convention of z = corresponding to the 
present epoch. 

Although the simpUcity of the above expression is attractive, 
we find that in some cases it provides a rather poor fit to the 
individual MAHs. We generalize the fitting formula by the fol- 
lowing two-parameter function 

M(fl) = 5^exp [<5 (1 - 1 /a)] , (4) 

The function in Eq. (2) is a special case of Eq. (4) for p = 0. 
The fits in which both a and p were varied and the fits with 
fixed p = are shown in Figure 1. These fits were obtained 
by minimization, even though the robustness of their rel- 
ative quality with respect to the choice of merit function was 
tested. For CL9 and 13 the fits with p = are rather poor. For 
the two-parameter fits to the MAHs of these clusters the value 
of a is close to zero, which means that the MAHs are better 
described by a power law in a, rather than by an exponential 
in I /a. Detailed study of the MAH shapes clearly requires a 
larger sample of halos. Analyzing the merger histories in terms 
of number of major mergers indicates that the power law behav- 
ior may be related to the high frequency of major mergers up 
to the present epoch. We note that the galaxy-size halos studied 
by W02 formed earlier on average than our halos, and thus only 
a small fraction (< 5%) of their halos were similar to our CL9 
and 13. 

Clearly, a smooth fitting function for the MAHs cannot cap- 
ture all the features of the actual evolution, such as minor and 
major merger events. These events however do influence the 
values of the best fit parameters. In the case of CL7, one can 
see the effect of a very recent major merger. As can be seen in 
Figure 1, if we make the fit at ao = 0.95 instead of ao = 1 (i.e., 
with the observation epoch prior to the merger), we get a better 
overall fit. 

van den Bosch (2002) found that the average mass accretion 
history of halos generated using the extended Press-Schechter 
(EPS) formahsm is described well by a different two-parameter 
function 

log(a) 
log(l+z/) 

with Zf, V the parameters to be determined. By definition, z/ is 
the redshift which corresponds to (M/Mq) = 1 /2. Typical best 
fit values for v are in the range of 1.4-2.3. As before, clusters 
CL9 and 13 are exceptions with u^l (i.e., their MAH is power- 
law M « M()a''). We find that Eq. (5) gives fits equally good to 
those obtained with the two-parameter function of Eq. (4) in all 
cases. This function however is not convenient when used to 
calculate the formation redshift via the logarithmic derivative 
of the mass with respect to the scale factor. More specifically, 
in Eq. (5) (and its derivative) z/ is by construction positive, and 
thus for a general i^, a has to be < 1. This will not be the case 
for objects whose mass accretion rate reaches 5 = 2 in the future 
(negative formation redshift), and thus Eq. (5) cannot be used 
in these cases to obtain a formation redshift. We choose to use 
Eq. (4), as it is a simple extension of the function used by W02 
and the definition of z/ via the mass logarithmic derivative is 
guaranteed. 



log(M)=-0.301(-l)^ 



(5) 
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Table 2 

Halo parameters from mass accretion histories and density fits 



Halo 


Zlmm 


IIM/M 


Zl/2 


2/ 


Best fit 


C-2 


slope 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


CLl 


1.24 


0.65 


1.09 


0.99 ±0.24 


M 


9.7 


-1.23±0.19 


CL2 


0.85 


0.56 


0.83 


1.37±0.14 


JS 


9.6 


-1.67±0.15 


CL3 


0.75 


0.36 


0.67 


0.78 ±0.19 


NFW 


10.7 


-1.35 ±0.20 


CL4 


0.95 


0.77 


0.78 


0.83±0.21 


JS 


5.4 


-1.89 ±0.20 


CL5 


0.85 


0.34 


0.69 


1.32±0.33 


M 


12.2 


-1.30 ±0.20 


CL6 


0.95 


0.45 


0.91 


1.33±0.33 


JS 


14.7 


-1.83±0.19 


CL7 


0.03 


0.82 


0.07 


1.04 ±0.25 


JS 


13.6 


-2.01 ±0.22 


CL8 


1.74 


1.83 


0.41 


0.80±0.20 


M 


11.5 


-1.25 ±0.22 


CL9 


0.15 


0.95 


0.19 


-0.04 ±0.14 


M 


3.5 


-1.68 ±0.22 


CLIO 


0.25 


0.78 


0.22 


-0.11 ±0.08 


M 


2.3 


-1.78 ±0.28 


CLll 


0.85 


2.11 


0.80 


0.21 ±0.09 


NFW 


8.4 


-1.38±0.31 


CL12 


1.24 


1.67 


1.23 


1.64 ±0.43 


M 


12.6 


-1.50 ±0.24 


CL13 


0.03 


0.33 


0.25 


0.41 ±0.16 


JS 


4.3 


-1.36±0.42 


CL14 


0.95 


0.81 


0.93 


0.82 ±0.20 


JS 


10.8 


-1.42±0.30 



Note. — (2): redshift of last major merger, (3): fractional mass change during last major merger, (4): redshift wliere iialf of tlie cluster's current mass has been 
accreted, (5): formation redshift defined from MAHs, (6): best fit among the Navarro et al. (1996, 1997), the Moore et al. (1998), and the Jing & Suto (2000) profiles 
(7): concentration index, (8): logarithmic slope as obtained by averaging the local logarithmic slope between the smallest resolved radius and 3% of the virial radius 



The formation redshifts estimated using the best fit param- 
eters of Eq. (3) and Eq. (4) are given in column 5 of Table 2. 
Note that this definition of Zj allows for future (negative) for- 
mation redshifts. The two definitions of the formation redshift, 
Zi/2 and z/, are correlated at 98% probability level (a Spearman 
rank correlation of 0.58). In addition, to evaluate the effect of 
the early and late portions of the MAH on the formation red- 
shift, we estimated the formation redshifts using only the parts 
of the MAH for which a < 0.65 and a > 0.65: zf'^'^ and z}'^'^^. 
The best fit values of both "'^^ and Zf^'^^ are consistent with 
the values of Zf within errors. In addition the values of z^"^^ 
are consistent with Zf^'^^ within (large) errors. Although our 
cluster sample is small, the significant spread in the values of 
Zf for clusters of the same Miso is apparent. This indicates that 
halos of the same mass exhibit a wide range of MAH shapes. 

5. density profiles 

5.1. The fitting procedure 

For each cluster halo, we fit the Navarro et al. (1997, NFW), 
the Moore et al. (1998, M), and the Jing & Suto (2000, JS) 
analytic density profiles. For a general profile of the form 

Ps 



(6) 



x'y(l+x")(/3-^)/"' 

the NFW, M, and JS profiles have values of {a, (3, 7): (1, 3, 1), 
(1.5, 3, 1.5), and (1, 3, 1.5), respectively. 

In what follows, we will define the concentration of a halo as 
c-2 = rm/f-i and Cv = r^lr-2, where r-2 is the radius where 
the logarithmic slope of the best fit profile is equal to -2, rjgo 
is the radius within which the average density is equal to 180 
times the mean matter density of the universe, and r^-^ is the 



virial radius defined using the redshift-dependent virial over- 
density (« 180 at z > 1 and « 340 at z = 0). To convert from 
(Ml 80, ri 80) to (Mvir, ^vir) wc usc the fitting formulas of Hu & 
Kravtsov (2003). For the general profile of Eq.( 6), the radius 
r-2 is given by 



r-i = 



7- 



2-/3 



l/a 



(7) 



where r^ is the scale radius of the corresponding analytic pro- 
file. Thus, r-2 = r, (NFW), r-2 « 0.63r, (M), r-2 = 0.5r, (JS). 
In other words, c_2 = cnfw, c-2 ~ 1 .59cm, and c_2 = 2cjs, if the 
best fit is found to be the NFW, the M, or the JS profile, respec- 
tively. We present the resulting concentration indices at z = in 
colunm 7 of Table 2. 

There is a number of factors that may affect the fits of ana- 
lytic profiles to the profiles of the simulated clusters: the choice 
of binning, the merit function, the range of radii used in the fit- 
ting, the weights assigned to the data points, etc. For example, 
for the merit functions sensitive to the number of bins, such as 
X^, the choice of binning and bin weights are extremely impor- 
tant. 

In the following analysis, we use equal-size logarithmic bins 
in order to give more statistical weight to the inner regions of 
the halos. The number of bins is thirty for late epochs. For 
early epochs the number of bins is reduced to ensure that each 
bin contains a sufficiently large number (> 100) of particles. 
We take as bin center the average radius of all particles in a bin. 
We checked that the fits are robust when varying the number 
of bins around the adopted value. We find, however, that the 
choice of binning affects the quality of the fit. For example, 
for a large number of bins the resulting profiles are quite noisy. 
Our choice of binning minimizes the noise. We weight the data 
points by the poisson noise in the number of particles of each 




Fig. 2. — Middle panel: NFW (dotted lines), Jing and Suto (short-dashed 
lines), and Moore et al. (dot-dashed lines) fits to the density distribution of 
CL2 (solid lines). The upper set of lines corresponds to results obtained by a 
minimization and the lower set (displaced by a factor of 10 for clarity) to 
fits obtained by minimizing the maximum absolute fractional deviation of the 
fits from the data. Upper panel: Fractional deviation of the fits (pf) from the 
data (pd) for fits obtained via minimization. Lower panel: Same as in the 
upper panel but for fits obtained via maximum absolute fractional deviation 
minimization. 



bin. 

The presence of substructure may substantially bias fits of 
smooth analytic profiles. In particular, substantial amount of 
substructure is present in the outer regions of halos and the 
profiles in these regions are often non-monotonic exhibiting 
"bumps". To minimize the bias, we fit the profiles using only 
the bins from a minimum resolved radius (see § 5.2) up to the 
radius within which the average density is equal to 500 times 
the critical density of the universe, rsoo- This choice is moti- 
vated by the results of Evrard et al. (1996) who find that the 
material within this radius is generally relaxed and in hydro- 
static equiUbrium. We find that for the clusters in our sample 
'■50o/''i80 — 0.36-0.37 at z = 0. To the same end, the density 
profiles of CL9 and 10 were obtained by averaging the z = 
and z « 0.05 outputs. This renders the profiles less noisy and 
improves the quahty of the obtained fits. The averaging does 
not change the best fit parameters significantly. 

It is important to understand that the formal quality of the 
fit may depend on the merit function, as well as the kind of 
binning and weighting used. In the following analysis, we fit 
the analytic profiles for the parameters ps and by minimiz- 
ing the x^- Klypin et al. (2001) show that merit function 
applied to the profiles with logarithmic binning with the Pois- 
son error weights results in the fits of the NFW profile that are 
systematically below the simulated profiles in the inner regions. 
The logarithmic biiming gives higher density of data points at 
small radii, creating thus bins with smaller number of parti- 
cles. The fits for the CL2 profile are shown in Figure 2. 



Also shown are the fits obtained when using the maximum frac- 
tional deviation (MFD) merit function, max[|pfit-pdata|//3data], 
which gives equal weight to all radial bins. Figure 2 shows that 
this merit function reduces the deviations in the inner region 
(r < 0.02ri8o) at the expense of significant deviations at inter- 
mediate radii (0.02 < r/rm < 0.2). Although the MFD merit 
function is less sensitive to the choice of biiming, it is more sen- 
sitive to the presence of substructure bumps in the profile than 
X^. Both merit functions have their pluses and minuses. 

Luckily, we find that regardless of the merit function and bin 
weighting used, the relative goodness of fits for different ana- 
lytic profiles remains the same. If, for example, the M profile is 
a better fit to the simulated profile than the NFW and JS profiles 
in the x^ minimization, it is the resulting best fit in the max- 
imum deviation minimization as well. The conclusion about 
which profile fits best is therefore robust. Note, however, that 
this is not true for the conclusions about the systematic ways by 
which a given fit fails. For example, the characteristic 'S' shape 
of the fractional deviation as a function of radius for the NFW 
fits found in various studies (e.g., Moore et al. 1999b; Ascasi- 
bar 2003) is not a robust result because it depends on the merit 
function, binning and weighting. 

In addition to the three widely used analytic profiles, we ex- 
perimented with fits of more general analytic expressions of 
the form given by Eq. (6). Overall, the fitting procedure with 
all three parameters a,/3,7 varied leads to strong degeneracies 
between parameters (see also Klypin et al. 2001). One can 
find several combinations of parameters that fit the data equally 
well. For example, good fits with inner asymptotic inner slopes 
as shallow as 7 = -0.3 can be found. 

In view of these degeneracies, we choose not to use general- 
ized fits but simply complement the analytic fits with measure- 
ments of the logarithmic slope profiles s{r) = d log p{r)/d log r. 
This analysis is complementary to the fits because the slope is 
sensitive to the local shape of the profile, while the fits may be 
sensitive to its global shape. The logarithmic slope is computed 
using the linear fit to log p- log r locally. We use five neigh- 
boring profile bins centered on a given bin in the fit (i.e., two 
bins on either side) with a total of 100 bins for the whole range 
of radii from r„,„ to rigo- The choice of the number of bins is 
a tradeoff between the slope errors and spatial resolution. We 
experimented with fitting polynomials up to the 4th order but 
found no advantage over a simple Unear fit. The local loga- 
rithmic slope is sensitive to the presence of transient massive 
substructures within the halo. For illustrative purposes only, to 
reduce the substructure-induced noise, we additionally smooth 
the slope using a tophat filter. 

5.2. Convergence Study 

To study the effects of mass and force resolution, CL2 was re- 
simulated with eight times more particles (mp = 3.95 X 10'/i"'Mq) 
and with more refinements. The ART code performs mesh re- 
finements when the number of particles in a mesh cell exceeds 
a specified threshold. Thus, the mass resolution is tightly linked 
to the peak spatial resolution achieved in simulation. The cell 
size of the highest refinement level, which we will consider to 
be the formal resolution of the simulation, was O.6/7"' kpc and 
1 .2/i~' kpc in the higher- (HR) and lower-resolution (LR) sim- 
ulation, respectively. The HR simulation was initiaUzed using 
the same set of modes as the LR. We therefore follow the forma- 
tion of the same object with more particles. The comparison of 
density profiles allows us to check for the two-body relaxation 
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Fig. 3. — Density profiles of the cluster CL2 in the low- (LR, solid line) 
and high-resolution (HR, dotted line) simulations {middle panel). The profiles 
were obtained by averaging the z = 0,0.02 and 0.1 outputs of the correspond- 
ing runs. The bottom panel shows the fractional deviation between the LR and 
HR profiles. The error bars are computed by propagating the shot noise in the 
density profiles. The top panel shows the local logarithmic slope as a func- 
tion of radius in the HR {squares) and the LR {triangles) runs. In the middle 
panel, the vertical arrow at ~ lOA"' kpc (or four times the formal resolution 
of the LR run) denotes the minimum distance used in our analyses. The ar- 
rows at large scales denote the radii corresponding to various commonly used 
overdensities: rsoo, the radius within which the average density equals 500 the 
critical density, and r34o and rigo, the radii within which the average density 
equals 340 and 180 times the mean density of the universe, respectively. 



effects, which may be important in cluster cores (Diemand et al. 
2003), and numerical convergence. 

We compare the density profiles of CL2 in the HR and LR 
simulations in Figure 3. In order to minimize the differences 
due to substructure, the profiles shown are obtained by aver- 
aging the z = 0, 0.02, and 0.1 outputs of the corresponding 
runs. The figure shows that the fractional difference between 
the profiles is < 0.2 down to ~ 3 formal resolutions of the LR 
run. This is in agreement with a previous convergence study 
for the ART code using simulations with lower mass resolution 
(Klypin et al. 2001). Comparison of density profiles of clus- 
ters in ART simulations with the density profiles in simulations 
using the Gadget code (Springel et al. 2001) was recently per- 
formed by Ascasibar et al. (2003), who found excellent agree- 
ment between the two codes at the resolved scales. 

The upper panel of Figure 3 shows the local logarithmic slope 
of the density profiles as a function of radius (see § 5.1 for de- 
tails). The error bars are computed by propagating the Poisson 
errors in the density profiles. At r > lOOh"^ kpc the strong non- 
monotonic variations of the slope are due to the presence of 
substructure. Despite the averaging, the small differences in the 
locations of substructures result in large differences in the slope 
value at a given r. At the same time, the slopes in the HR and 
LR runs agree well at scales 5 < r < IQQhr^ kpc. It is interesting 
to note that there is no evidence for a well-defined asymptotic 




Fig. 4. — The density profile of CL2 at (from top to bottom) 7. = 0, 0.2, 0.4, 
1, and 1 .5 {solid lines), for the low (LR, left panel) and high (HR, right panel) 
resolution runs. The profiles at z > are scaled down by a factor of 10 with 
respect to each other Also shown are the best fit NFW {dotted lines) and the 
JS {dashed lines) profiles. The figure shows that at all shown redshifts the JS 
profile is a better fit to the simulated profiles than the NFW in both runs. 



inner slope. The local logarithmic slope in both runs mono- 
tonically increases with decreasing radius down to the smallest 
resolved scales. 

Based on these results, in the subsequent analysis we will 
conservatively consider only scales greater than rnu„ = \Qh~^ co- 
moving kpc or (eight formal resolutions of the LR run) for both 
the fits and the plots. In addition, we require that more than 
200 particles are contained within the minimum radius (Klypin 
et al. 2001). In our simulations this criterion is relevant only at 
early epochs (z > 2), since at later epochs a \0h~^ kpc radius 
always contains more than 200 particles for all clusters. 

One of the main results of our analysis is an apparent diver- 
sity of the density profiles. At the same time, we find that if a 
profile is best fit by a particular analytic profile at z = 0, it is gen- 
erally best fit by the same analytic profile at early epochs out to 
z ^ 1-2. Potentially, this is an interesting clue to the processes 
that determine the shape of the profile. It is therefore important 
to check that the conclusion does not change with resolution. 
Figure 4 shows the fits of the NFW and JS profiles to the pro- 
files of CL2 in the LR and HR runs at different epochs. The fits 
were done using bins in the radial range [rmin,r5oo] (see § 5.1). 
The figure shows that at all shown redshifts the JS profile fits 
the simulated profiles at small radii better than the NFW pro- 
file in both runs. This is remarkable as the cluster experiences 
fairly rapid increase in mass and several violent mergers be- 
tween z = 1 .5 and z = 0. The mass changes by more than a 
factor of five during this period (see Fig. 1). The cluster un- 
dergoes an intense merger event at z ~ 0.6 and the last major 
merger for our definition occurred at z = Zlmm — 0.85. 

We find that our fitting results are robust to changes in both 
the minimum and maximum radius used in the fits. For exam- 
ple, concentration changes by no more than ~ 10 — 20% if a 
different outer radius is used (~ 2/3ri8o, and for some clusters 
an outer radius ~ riso did not change the results much). More 
importantly, the conclusion about the best fit analytic profile 
remains the same, although in some cases the best fit changes 
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from the M to the JS or vice versa. This can be expected be- 
cause these analytic profiles are quite similar We also repeated 
fits with twice as large minimum radius (« 20h~^ kpc comov- 
ing) and find that the best fit analytic profile remains the same 
and that the concentration changes by < 10%. 

5.3. Results 

Our results on the density profiles at z = are presented in 
Figure 5. In the middle panel we plot the actual profiles as well 
as the NFW, M, and JS fits. All profiles are plotted atr < r^^o 
but are fit using only bins in the range rmin < r < rsoo, as dis- 
cussed in § 5.2. In the lower panel we present the fractional 
deviation as a function of distance for each of the analytic fits. 
The analytic profile that provides the best fit to the profile of 
each cluster is given in column 6 of Table 2. In the upper panel 
we plot the logarithmic slope as a function of radius for each of 
the analytic fits and the actual local logarithmic slope as calcu- 
lated from the simulated profiles (see § 5.1). 

In most cases the best analytic fit provides by far the best fit 
compared to the other profiles. This is especially true for CL3 
and 1 1 (see Fig. 5), the two clusters best fit by the NFW profile. 
The other two analytic profiles fail significantly compared with 
the NFW for these clusters. If we consider the similar M and 
JS as one family of profiles, the two families (NFW and M/JS) 
typically differ significantly in quality. As discussed above, the 
systematic way in which the NFW profile fails to fit the data can 
be attributed to the merit function used to obtain the fits. For our 
choice of merit function, the fits follow the actual profile well 
at intermediate distances. The largest deviations occur at the 
innermost regions. For large distances, the three fits are almost 
indistinguishable. 

Figure 5 and Table 2 clearly show that there is significant dis- 
persion in the shapes of the profiles, concentrations, and inner 
slopes. The dispersion of concentration parameter was stud- 
ied in several analyses (Navarro et al. 1997; Jing 2000; Bullock 
et al. 2001; Eke et al. 2001; Wechsler et al. 2002; Zhao et al. 
2003b) and is thought to be related to the distribution of the 
halo formation epochs (Wechsler et al. 2002). The typical val- 
ues of scatter are fTioge^ ~ 0.14 with only a weak dependence 
on mass(Bullock et al.' 2001; Wechsler et al. 2002). Table 2 
shows that the c-2 concentration indices of our clusters span a 
wide range of values, from 2.3 to 14.7. Making the appropriate 
conversion from c-2 to Cj,, we find that the formal dispersion is 
ciogcv ~ 0-2 at z = 0, larger than the dispersion for the smaller 
mass halos used to derive this scattering in other studies, which 
may reflect the more recent formation times of cluster halos, as 
well as their more diverse MAHs. For example, Klypin et al. 
(2003) and Colin et al. (2003) find a significantly smaller dis- 
persion cTioge, « 0.1 for a subsample of relaxed halos without 
significant substructure. Indeed, the clusters with the three low- 
est concentrations, CL9, 10, and 13 have all had a recent major 
merger (see zlmm in column 2 of Table 2). The small concen- 
tration of these objects is due to the shape of their density pro- 
files which are close to a power law over a wide range of radii. 
Careful examination of MAHs and merger histories, indicates 
that the recent major merger activity results in a low concentra- 
tion of density profiles. CL4, which also has a small concen- 
tration compared to the majority of the clusters, has a formal 
Zlmm = 0.95, for the definition of major merger adopted in our 
study. However, in agreement with the other low concentration 
objects, it had a large merger (~ 22% fractional mass increase) 
at a very recent epoch {z ^ 0.15). In addition, we find a strong 



correlation between z/ and c-2, in agreement with the correla- 
tion advocated by W02. 

Figure 6 shows the redshift evolution of the median virial 
concentration, Ci,, of our sample. We also plot the predictions of 
the models by Bullock et al. (2001) and Eke et al. (2001). Our 
results seem to be in some general agreement with both model 
predictions. Overall, the Eke et al. (2001) seems to be in bet- 
ter agreement than the Bullock et al. (2001) model. Recently, 
Dolag et al. (2003) found fliat the standard Bullock et al. (2001) 
recipe systematically underestimates concentrations of cluster- 
size halos at all redshifts. Figure 6 shows a similar trend in 
our simulations, although the difference we find is noticeably 
smaller This may be due to the smaller mean mass of clus- 
ters in our sample. Zhao et al. (2003b), for example, show that 
discrepancy between simulations and the BuUock et al. (2001) 
model increases with increasing halo mass. In the mass range 
probed here the difference from the analytic prescription of Bul- 
lock et al. (2001) is considerably smaUer than the scatter in con- 
centrations. 

Figure 7 shows the evolution of the average concentration 
and the average MAH of our clusters. The figure shows that 
during the period of rapid mass growth the concentration is 
nearly constant at Cy ~ 3 - 4, while during the period of grad- 
ual mass growth it increases with decreasing redshift as Cy cx 
(1 +zT^ . Therefore, the concentration of halos is approximately 
constant while they experience frequent maj or mergers and there 
may exist a "floor" to the concentration values of Cmin ~ 3, 
while concentrations start to increase with time at z < Zf . This 
behavior was pointed out by Zhao et al. (2003b,a). However, 
unlike Zhao et al. (2003b), we find that the model of Wechsler 
et al. (2002) does not underestimate the concentrations at the 
cluster-size halos of our sample. The evolution of average con- 
centration in Figure 7 is similar to that found by Dolag et al. 
(2003). The mean concentration of cluster progenitors in their 
simulations (their Figures 4 and 7) is approximately constant at 

We do not find a clear connection between the redshift of last 
major merger and the best fit profile (column 6 in Table 2). We 
do find that none of the objects with low zlmm has the NFW pro- 
file as the best fit, but the statistics of our sample is too small to 
reach a firm conclusion. Nevertheless, for each individual sys- 
tem the shape of its density profile is remarkably stable during 
evolution. The best fit analytic profile at z = is typically also 
the best fit at earher epochs, as was shown for CL2 in Figure 4. 
Figure 8 shows evolution of the density profiles and the analytic 
fits at each epoch for CLl and CL3. The NFW profile is a better 
fit than the JS at all epochs for z < 1-5 for CLl. This stability 
of the profile shape with time holds for most clusters with some 
exceptions. CL3 illustrates the case where the best fit analytic 
profile changes from epoch to epoch. We find this behavior for 
four out of the fourteen clusters in our sample. 

Note that the NFW fits never reach their inner asymptotic 
slope of -1 at the radii we probe. The average logarithmic 
slopes estimated by averaging the local slope around 0.03rigo 
(see §5.1) range from -1 .2 to -2, and are given in column 8 of 
Table 2. In most cases, the local slope changes monotonically 
with radius with no sign of reaching the asymptotic iimer slope. 
Note that for typical concentrations of cluster halos we expect 
the asymptotic slope to be reached at the resolved scales, at least 
for the M profile^. This can be seen from the slope profiles for 

' From Eq. (6) is the radius where the logarithmic slope is equal to -(/3+7)/2 
with the asymptotic slope reached at r « r^. The NI'W and the JS profiles 
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Fig. 6. — Median concentration vs. virial mass at different redsliifts for tlie 
progenitors of clusters in our sample (points). The vertical errorbars represent 
the la scatter in concentrations for the 14 clusters, while horizontal errorbars 
show the mass range of the halos at each epoch. The predictions of the Bullock 
et al. (2001) (thick lines) and the Eke et al. (2001) (thin lines) models are 
plotted for comparison. 




Fig. 7. — Average mass accretion history (top panel) and average concentra- 
tion of cluster progenitors (bottom panel) as a fimction of scale factor measured 
in units of the formation scale factor, at. The errorbars in both panels repre- 
sent the 1<T spread aroimd the mean. The figure shows that the concentration 
is approximately constant at Cv 3 - 4 during the period of rapid mass ac- 
cretion (a /at < 1) and increases with decreasing redshift during the period of 
slow mass growth (a/at > !)• For comparison the dashed line shows (1 +z)~' 
evolution. 



have typically smaller rj than that of the M profile. As a result, they reach their 
asymptotic slopes at smaller distances. An analytic profile, of course, can be a 




Fig. 8. — The density profile of CLl (left panel) and CL3 (right panel) 
at (from top to bottom) z = 0, 0.2, 0.4, 1, and 1.5 (solid lines). The profiles 
at z > are scaled down by a factor of 10 with respect to each other. Also 
shown are the best fit NFW (dotted lines) and the JS (dashed lines) profiles. 
At all shown redshifts the NFW profile is a better fit to the simulated profile of 
CLl than the JS profile. CL3 is shown as an exception to the typical behavior 
represented by CLl and CL2 (Fig. 4). For CL3 the best fit at z = (NFW) is 
not the best fit at earlier epochs. 



the best fit analytic fits shown in the top panels of Figure 5. 
The only cases where the slope profiles are flat for relatively 
large radius ranges, correspond to the halos with recent merger 
activity and rapid mass growth (i.e., CL4, CL9, CLIO, CL13). 

Interestingly, we find that the density profiles of systems that 
experience intense merger activity until the present epoch (clus- 
ters CL4, 9, 10 and 13) can be well described by a single power 
law r"^ with slope 7 ranging from -1.5 to -2. Similarly to other 
profile shapes, the power law density profile for these systems 
is maintained for earUer epochs out to z ~ 1.5. In addition, 
there is evidence that the profiles of all clusters during their 
rapid mass growth stages are close to a power law. In particu- 
lar, we find that the power law provides an increasingly better fit 
with increasing redshifts for all of our clusters. At early epochs 
(z > 1.5-2), the power law fit is always either comparable or 
better than the NFW, M, and JS analytic profiles. The power 
law like profiles relate to low concentrations. This is because 
they maintain a slope slightly shallower than -2 out to large 
radii so that the scale radius is large. Indeed, the clusters with 
power law profiles at z = have the lowest concentrations in 
our sample. The decrease of concentrations at higher redshifts 
may thus reflect the power law like density profiles of actively 
merging systems. 

6. DISCUSSION AND CONCLUSIONS 

We have studied the MAHs and density proflles of 14 cluster- 
size halos simulated using the ART code in a flat ACDM cos- 
mology. In agreement with previous studies, we find that most 
MAHs have a similar shape: an early, merger-dominated mass 
increase followed by a more gradual, accretion-dominated growth. 
To obtain a formation redshift that characterizes the overall shape 
of the MAH we perform analytic functional fits. The typical 
MAHs are well described by the one parameter exponential 

good fit regardless of whether its iimer asymptotic slope is resolved. 
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function proposed by W02 (Eq. [2]). Two of the clusters in 
our sample experience intense merger activity and rapid mass 
growth until the present-day epoch. The MAHs of these sys- 
tems are better described by a one parameter power-law func- 
tion in the scale factor. We thus generalize the form proposed 
by W02 into a two parameter form (Eq. [4]) to encompass both 
exponential and power-law MAHs. For each class, however, 
the fit reduces to a one parameter fit. 

We check the convergence of halo density profiles using a re- 
simulation of one of the clusters with eight times more particles 
and better force resolution. We show that both the halo profiles 
and their local logarithmic slopes converge at scales larger than 
about four times the formal resolution of the low resolution run, 
in agreement with a previous convergence study of the ART 
code by Klypin et al. (2001). We fit the density distribution of 
the clusters with the NEW, M, and JS analytic profiles. Exper- 
iments show that the choice of merit function, weighting, and 
binning affect the absolute quality of a fit and may bias conclu- 
sions about how well a particular analytic profile fits simulation 
results. We find, however, that the relative goodness of fit for 
the three analytic profiles and our conclusions about the best fit 
profile are robust to the changes in binning and merit function. 

The main general result of our study is a remarkable diver- 
sity of the mass accretion histories, profile shapes, concentra- 
tions, and iimer slopes for cluster-size halos in a relatively nar- 
row mass range. The concentrations of cluster-size halos at the 
present-day epoch exhibit a scatter of crioge,, ~ 0.20. This scat- 
ter is related to the diversity of halo MAHs and formation red- 
shifts. We find a statistically significant correlation between 
the formation redshift and the concentration of a halo, in agree- 
ment with results of W02. There is a more detailed connec- 
tion between the MAH and concentration. The concentration 
of a halo is approximately constant at Cv w 3 - 4 during the pe- 
riod of rapid mass growth and frequent major mergers (z > Zf) 
and increases with decreasing redshift when the mass accre- 
tion rates slows down at z < Zf- This behavior was recently 
pointed out by Zhao et al. (2003b, a). The implied "floor" in the 
concentration is not accounted for in the currently used models 
for Cv(M) (BuUock et al. 2001; Eke et al. 2001), which predict 
a monotonic decrease of concentration with increasing mass 
and may thus underestimate concentrations of the most mas- 
sive, > 5 X lO'^^/i"^ Mq, halos. This may have important im- 
plications for estimates of expected number of wide- separation 
quasar lenses (e.g., Kuhlen et al. 2003) and other results sensi- 
tive to the concentrations of very massive clusters. 

The inner logarithmic slope of cluster profiles at 3% of the 
virial radius (or 10-50 kpc) ranges from -1.2 to -2. In the 
best resolved clusters the logarithmic slope does not seem to 
reach a specific asymptotic value down to the smallest resolved 
scales in our simulations (r/ri8o ~ 0.007). A similar conclu- 
sion was reached by Klypin et al. (2001) for galaxy-size halos 
and several recent studies (Power et al. 2003; Ascasibar et al. 
2003; Fukushige et al. 2003; Hoeft et al. 2004; Hayashi et al. 
2003). It is still not clear whether the density profiles in our 
simulations are consistent with density distribution of observed 
clusters. We note, however, that at the scales probed in obser- 
vations the slope is not expected to be shallower than -1. 

The asymptotic value of the slope has been a subject of much 
numerical effort in the last several years. Our results indicate 
that a universal asymptotic slope may not exist. We should 
note that the resolution of current dissipationless simulations 
is sufficiently high to converge on the density profile at scales 



smaller than the size of a typical central galaxy in clusters and 
groups (~ 30-50 kpc). Further improvement in profile mod- 
eling should therefore include realistic dynamics and cooling 
of the baryonic component as contraction of gas is expected to 
significantly modify dark matter distribution at these scales. 

One of the most interesting results of our study is existence 
of systems with density profiles that can be well described by 
a power law p oc with 7 ranging from « -1 .5 to « -2. All 
of these systems are still in their rapid mass growth stage and 
experienced a recent major or minor merger. Remarkably, these 
halos maintain the power law density profiles at earUer epochs 
out to at least z ~ 1 .5. The relatively shallow 7 > -2 power law 
slopes result in low concentrations as the scale radius where the 
density profiles reaches the slope of -2 is at large radii. There 
are also indications that the profiles of all clusters are power 
law Uke during their rapid mass growth stages. We find, for 
example, that the power law provides an increasingly better fit 
with increasing redshifts for all of our clusters. At early epochs 
(z > 1.5- 2), the power law fit is always either comparable or 
better than the NEW, M, and JS analytic profiles. We did not 
find any correlation of the power slope with the details of the 
cluster MAH. It would be interesting to look for such correla- 
tions using a larger sample of objects. 

When the mass growth slows down at z > Zf , an outer steeper 
density profile is built up. As pointed out by Zhao et al. (2003b), 
the difference in density profiles during the two mass accre- 
tion regimes may be due to more violent and thorough relax- 
ation during the period of rapid mass growth. Although Zhao 
et al. (2003b) focused on halo concentrations and did not con- 
sider density profile shapes, their results are consistent with our 
conclusions. In particular, they find that during the rapid mass 
growth stage the circular velocity is nearly constant from the 
scale radius to the virial radius. This behavior is consistent with 
a power law density distribution with a slope close to -2. 

In a recent study, Ascasibar et al. (2003) find that objects 
which experienced a recent merger event* have lower concen- 
trations and steeper inner profiles than more relaxed systems. 
This is consistent with our findings described above. At the 
same time, Ascasibar et al. (2003) find that relaxed systems 
are better fit by the NEW, while systems with a recent major 
or minor merger by the M profile (see also Ascasibar 2003). 
They thus associate a particular shape of the profile with a re- 
cent merger history. In contrast, our results show the shape of 
density profiles is set early in the halo evolution and is usually 
stable over the past ten billion years. Clusters with density pro- 
files best described by the NEW rather than a JS at z = 0, tend 
to have NEW-like profiles at earlier epochs as well. The reverse 
is also true. We tested this conclusion using the high-resolution 
re-simulation of one of the clusters in our sample. Also, we 
do not find any correlation between the redshift of last major 
merger (in our definition) or the formation redshift and the best 
fit analytic profile. 

The origin of the distinctive density profile shape of the CDM 
halos remains poorly understood. Our results and results of 
other recent studies indicate that the shape is tightly linked to 
the halo mass accretion history. During the period of rapid 
mass accretion the violent relaxation is significant and results 
in a power-law Uke density distribution. This stage of evolu- 

* Note that Ascasibar et al. (2003) identify a recent merger by the presence 
of massive substructures within virial radii of their systems. This is different 
from our definition, which identifies major mergers directly from mass accre- 
tion tracks. 
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tion usually occurs early when the universe is dense and builds 
up the inner dense regions of halo. At this point, the logarith- 
mic slope of the density distribution is shallower than -2 over a 
large fraction of the halo volume and its concentration is small. 
At later epochs, as the mass accretion rate slows down, the outer 
regions of the halo are built while its central regions remains 
nearly intact. This can be seen in Figure 8 and Figures 10-13 
of Fukushige & Makino (2001). This picture can explain why 
the best fit analytic profiles tend to be the same at various red- 
shifts. The fits are sensitive to the density distribution at small 
and intermediate (~ r^) radii which are set early. However, it is 
still unclear which process(es) determines a particular shape of 
the profiles. The key to understanding these processes appears 
to be in the details of early evolutionary stages of CDM halos, 
which will be the subject of a future study. 
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Fig. 5. — Middle panel: fits of the NFW (dotted curves), the Jing and Suto (short-dashed curves), and the Moore et al. (dot-dashed curx'es) profiles to the density 
distribution of the clusters of our sample at z = 0. The fits were done using the range [r,„in. ''soo] ™d a merit function (see §5.1 for details). The way the choice 
of merit function changes the fits can be seen in Figure 2. Bottom panels: deviations of each one of the fits (p/) from the actual profile (pd)- Top panels: local 
logarithmic slope as a function of radius for the 3 fits. The points correspond to the local slope as derived from the actual profile. 
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